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^ ' ■ Abstract 

I The paper is concerned with the problem of shape preserving interpolatory 

subdivision. For arbitrarily spaced, planar input data an efficient non-linear 
subdivision algorithm is presented that results in limit curves, reproduces 
conic sections and respects the convexity properties of the initial data. Signif- 
icant numerical examples illustrate the effectiveness of the proposed method. 
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^ 1 Introduction and state of the art 

Subdivision schemes constitute a powerful alternative for the design of curves and 
^ ' surfaces over the widely studied parametric and implicit forms. In fact, they offer 

• a really versatile tool that is, at the same time, very intuitive and easy to use and 

implement. This is due to the fact that subdivision schemes are defined via itera- 
tive algorithms which exploit simple refinement rules to generate denser and denser 
point sequences that, under appropriate hypotheses, converge to a continuous, and 
potentially smooth, function. 

In the univariate case, the iteration starts with a sequence of points denoted 
by p° = {Pi '■ « G Z), attached to the integer grid, and then for any k > one 
subsequently computes a sequence p'^"^"^ = Sp^, where S : i{Z) — )■ ^(Z) identifies 
the so-called subdivision operator. 

Subdivision operators can be broadly classified into two main categories: inter- 
polating and approximating [121 126]. Interpolating schemes are required to generate 
limit curves passing through all the vertices of the given polyline p°. Thus they 
are featured by refinement rules maintaining the points generated at each step of 
the recursion in all the successive iterations. Approximating schemes, instead, are 
not required to match the original position of vertices on the assigned polyline and 
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thus they adjust their positions aiming at very smooth and visually pleasing limit 
shapes. As a consequence, while in the case of approximating subdivision the refine- 
ment rules rely on a recursive corner cutting process applied to the starting polygon 
p°, in the case of interpolatory subdivision, in every iteration a finer data set p'^"'"^ 
is obtained by taking the old data values p*^ and inserting new points in between 
them. Every such new point is calculated using a finite number of existing, usually 
neighboring points. In particular, if the computation of the new points is carried 
out through a linear combination of the existing points, the scheme is said to be 
linear, otherwise non-linear. 

Then, inside the above identified categories, the schemes can also be further 
classified. More specifically, they can be distinguished between stationary (when the 
refinement rules do not depend on the recursion level) and non-stationary; between 
uniform (when the refinement rules do not vary from point to point) and non- 
uniform; between binary (when the number of points is doubled at each iteration) 
and N-ary, namely of arity N > 2. 

Most of the univariate subdivision schemes studied in the literature are binary, 
uniform, stationary and linear. These characteristics, in fact, ease to study the 
mathematical properties of the limit curve, but seriously limit the applications of 
the scheme. Exceptions from a binary, or a uniform, or a stationary, or a linear 
approach, have recently appeared (see for example [21] and references therein), but 
none of the proposed methods provides an interpolatory algorithm that can fulfill the 
list of all fundamental features considered essential in applications. These features 
can briefiy be summarized as: 

(i) generating a visually-pleasing limit curve which faithfully mimics the be- 
haviour of the underlying polyline without creating unwanted oscillations; 

(ii) preserving the shape, i.e., the convexity properties of the given data; 

(iii) identifying geometric primitives like circles and more generally conic sections, 
the starting polyline had been sampled from, and reproducing them. 

Requirement (i) derives from the fact that, despite interpolating schemes are consid- 
ered very well-suited for handling practical models to meet industrial needs (due to 
their evident link with the initial configuration of points representing the object to 
be designed), compared to their approximating counterparts, they are more difficult 
to control and tend to produce bulges and unwanted folds when the initial data 
are not uniformly spaced. Recently this problem has been addressed by using non- 
uniform refinement rules [6] opportunely designed to take into account the irregular 
distribution of the data. But, although their established merit of providing visually 
pleasing results, there is no guarantee that such methods are convexity-preserving, 
i.e., that if a convex data set is given, a convex interpolating curve can be obtained. 
This is due to the fact that, such non- uniform schemes are linear and, as it is well- 
known [17], linear refinement operators that are cannot preserve convexity in 
general. 

The property (ii) of convexity preservation is of great practical importance in 
modelling curves and surfaces tailored to industrial design (e.g. related to car. 
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aeroplane or ship modelling where convexity is imposed by technical and physical 
conditions as well as by aesthetic requirements). In fact, if shape information as 
convexity is not enforced, interpolatory curves, though smooth, may not be satisfac- 
tory as they may contain redundant wiggles and bumps rather than those suggested 
by the data points, i.e., they feature unacceptable visual artifacts. Preserving con- 
vexity, while a curve is interpolated to a given data set, is far from trivial. But 
much progress has been made in this field, evidence of which is given by the recent 
burgeoning literature. In most publications, the introduction of subdivision schemes 
fulfilling requirement (ii) has been achieved through the definition of non-linear re- 
finement rules. In fact, despite linear subdivision schemes turn out to be simple 
to implement, easy to analyze and affine invariant, they have many difficulties to 
control the shape of the limit curve and avoid artifacts and undesired inflexions that 
usually occur when the starting polygon p° is made of highly non-uniform edges. 
Non-linear schemes, instead, offer effective algorithms to be used in shape-preserving 
data interpolation P EB [131 [IZl [ISl EH] • 

On the basis of the well-known, linear Dubuc-Deslauriers interpolatory 4-point 
scheme [10] for example, several non-linear analogues have been presented in order 
to solve at least one of the three above listed properties. On the one hand, non- 
linear modiflcations of the classical 4-point scheme have been introduced to reduce 
the oscillations that usually occur in the limit curve when applying the reflnement 
algorithm to polylines with short and long adjacent edges. These have been pre- 
sented in [8] and [H], and as concerns the case of convexity-preserving strategies 
(which are the ones capable of completely eliminating the artifacts arising during 
the subdivision process), we flnd the papers [15] and [19]. On the other hand, for the 
purpose of enriching the Dubuc-Deslauriers 4-point scheme with the property (iii) of 
geometric primitives preservation, a non-linear 4-point scheme reproducing circles 
and reducing curvature variation for data off the circle, has been deflned [25]. With 
the same intent, another modiflcation of the classical 4-point scheme in a non-linear 
fashion, had been given in [3]. 

With these papers, the theoretical investigation of non-linear interpolatory sub- 
division has only begun. A lot is still to be done, in particular as concerns the 
use of non-linear rules for reproducing salient curves other than circles, considered 
of fundamental importance in several applications. So far, it has been shown that 
non-linear updating formulas can be used in the deflnition of non- stationary subdi- 
vision schemes aimed at reproducing polynomials and some common transcendental 
functions. In particular, [23] respectively [U El [22] present subdivision algorithms 
that turn out to be circle-preserving respectively able to exactly represent any conic 
section. While the flrst is able to guarantee reproduction starting from given sam- 
ples with any arbitrary spacing, for the latter ones the property of conies precision 
is conflned to the case of equally-spaced samples. Most recently a shape and circle 
preserving scheme has been presented in [H]. 

Therefore, an outstanding issue that should be considered is the possibility 
of deflning an interpolatory subdivision scheme that is at the same time shape- 
preserving and artifact free, as well as capable of generating conic sections starting 
from any arbitrarily-spaced samples coming from a conic. This is exactly the pur- 
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pose of this paper. Based on an approximation order four strategy presented in [T] 
for estimating tangents to planar convex data sequences, we propose a convexity- 
preserving interpolatory subdivision scheme with conic precision. This turns out to 
be a new kind of non-hnear and geometry-driven subdivision method for curve in- 
terpolation. In Section [2] we start by describing the refinement strategy which relies 
on a classical cross-ratio property for conic sections and uses the tangent estimator 
from [1], for the case of globally convex data. In Section [3] we adapt the scheme to 
general, not necessarily convex data by segmenting the given polygon into convex 
segments and by presenting new refinement rules next to the junction points. In 
Section H] we summarize the whole subdivision algorithm in all its steps. Then, in 
Section |5] we present an adaptive variant of the scheme which is aimed at producing 
regularly spaced points in every round of subdivision. Section |6] contains proofs for 
the scheme's shape preservation and conic reproducing properties, as well as for the 
continuity of its limit curve. Section [7] is devoted to illustrating the scheme by 
several significant apphcation examples, and we conclude in Section |S1 



2 Definition of the scheme for globally convex 
data 

In this section we define a convexity preserving subdivision scheme for globally 
convex data which will then be the basis for the final shape preserving subdivision 
algorithm for general data. 

Curve subdivision schemes iteratively apply a subdivision operator S* to a start- 
ing point sequence p" = (p° : i G Z) yielding a new sequence p'"'^^ = Sp^ for any 
level k > 0. Our scheme being interpolatory, has refinement rules of the following 
form: 

fc+i _ k 

P2t+1 = ^{Pi-u, Pt Pi+1, P'^+u, Phu+l, P) 

where z/ = 2 is the number of points taken into account in the left and right hand 
neighborhoods of the segment pfp^i^j^ in order to define the newly inserted vertex 
P2i+i) p is a parameter point specified later. y9 is a non hnear function, which 
we will define in form of an algorithm. 

In order to detail the idea of this convexity-preserving scheme, we thus consider 
the following problem where, for simplicity, we omit the upper indices k. 

Problem 1. Given n points Pi((pi)a;, {Pi)y), i = 1, ■ ■ ■ , n {n > 5) , in convex position 
in the affine plane, we wish to obtain one new point Uj related to the i-th edge pi Pi+i ■ 

We will carry out the construction of the new points in the projectively extended 
affine plane. To this end we denote the projective counterparts of the affine points 

Pi{{Pi)x,{Pi)y), i = by Pi(Pi,o,Pi,i,Pi,2) where pi^ = = {pi)x,Pi,2 = 

\Pi)y 
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By projective geometry's principle of duality, a line L may be represented either 
by a linear equation 

IqXq + hxi + I2X2 = 

in variable point coordinates {xq, Xi,X2) or by a triple {Iq, li, I2) of constant line coor- 
dinates. The line coordinates of the line L(/o,^i,^2) joining two points Xi(xi,o, 
Xi^2) and ^2(2^2,0) 3:2,1, 2:2,2) may simply be calculated by the vector product Xi A 
X2 = L. In the same way, the point coordinates of the intersection point P{po, P11P2) 
of two lines Liili^, /i^i, li^2) and L2{l2,o, h,i, ^2,2) is obtained as L1AL2 = P. Without 
loss of generality we apply a normalization to the homogeneous point coordinates 
such that xq e {0, 1} for all calculated points. 

In order to preserve global convexity of the points we apply the following pre- 
processing procedure to the point set. In every given point Pi we estimate a tangent 
from a subset of five points (including the point Pi) by the conic tangent estimator 
presented in [1]. If the given points represent a closed polygon, then the five-point 
subset is composed of the point Pi, its preceding two points Pj_i, Pi-2 as well as its 
successive two points Pi+i,Pj_|_2 (by considering Pq = Pn, P-\ = Pn-i, Pn+i = Pi, 
Pn+2 = ^2)- If the given points represent an open polygon, then for z = 3, . . . , n — 2 
we proceed as above and for i G {1,2} (respectively i G {n — l,n}) the points 
{Pi,P2,P3,P4,P5} (respectively {Pn-4, Pn-s, Pn~2, Pn~i, Pn}) are taken. 




Figure 1: Illustration of the tangent construction. 

In order to apply the conic tangent estimator from [1] we locally rename the five 
points around Pj by Q3 = Pj, and by arbitrarily mapping the remaining four points 
to Qi, Q2, Qi, Q5 by a one-to-one map. The desired tangent in the point Q3 is then 
calculated by the formula (see pQ) 

M33:=g3A(Mi5A(AAP)), (2) 

where Mij = Qi A Qj for i 7^ j, A = A M34, and B = M^^ A M32. See Figure [D 
for an illustration. 

We then denote the obtained line in the point Pj by Lj and intersect every two 
consecutive lines generating the intersection points 

Ti = Li A Lij^i , 
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Figure 2: Convex delimiting polygon (dashed line) for a given convex data set (solid 
line). 

see Figure [21 

The dashed lines in Figure [2] constitute a convex delimiting polygon for the new 
points generated in the next subdivision level. If the initial points come from a conic 
section, the constructed lines Lj are the tangents to this conic in the respective points 
Pi. Otherwise the lines Lj approximate the tangents with approximation order 4, 
see [Ij. 

After this preprocessing step we now get back to the initial subdivision problem 
[H i.e., between every two points Pi and Pj+i insert a new point Ui by applying 
a classical result from projective geometry which, for the readers convenience, we 
recall in the following theorem (see e.g., [21 [71 [21]). 

Theorem 2. a) Let X, E, Eq, Ei be four points of a projective line Vi, where 
the points E, Eo, Ei are mutually distinct, and let [xq, xi) be the projective 
coordinates of the point X with respect to the projective coordinate system 
{Eq,Ei;E} of Vi. Then, the cross ratio cr{X, E, Eq, Ei) of the four points 
X, E, Eq, El in this order is defined by 

cr{X,E,Eo,Ei) = ^. 

Xq 

b) Let Pi,P2 be two points on a conic section c, and ti,t2 the tangents of c in 
Pi, P2 respectively, and let T be the intersection point ofti and ^2 {T = ti ("1^2)- 
Then, the point T and the line P1P2 are pole and polar with respect to the conic 
c. Let's further denote the intersection points of any line It through T with 
the conic c by P and U , and the intersection point of It and T 's polar P1P2 
byX{lTnc= {P, ?7}, /t n P1P2 = X) . Then 

cr{U,P,X,T) = -l, 

and the four points {U, P, X, T) are said to be in harmonic position. 

Theorem [21 gives us the means of constructing a point U from three known 
collinear points P, X, T such that the harmonic cross ratio condition for conic sec- 
tions is satisfied. For an illustration see Figure [3l 



6 



Figure 3: Harmonic cross ratio condition for conic sections. 



For the choice of the parameter point P which is needed for the construction 
of the new point Ui inside the triangle APjTjPj+i, the whole region bounded by 
the lines Lj, Lj+i and PiPi^i containing the given convex polygon (see Figure S]) is 
suitable. In particular, every point Pj^ of the given convex polygon can be taken as 
parameter point (with exception of Pi and Pj+i), and since such a choice guarantees 
the reproduction of conic sections we opt for it. 




Figure 4: Suitable region for the choice of the parameter point P. 



Let's thus denote by Xj the intersection point of the lines PiPi+i and Pj^Ti for a 
chosen ji e {1, . . . ,n} \ {i,i + 1}, i.e., Xi = PiPi+i fl P^-.Tj. 

In order to guarantee a regular distribution of the inserted points we propose 
to choose the index ji by the following angle criterion. To this end we temporarily 
come back to the Euclidean plane and introduce the midpoint rrij for each segment 
PiPi+i. Let gi = tjiiij, hij = tipj be the connecting lines of the points tj and nij 
respectively tj and pj, and a* = Z{gi,hij) the angle between these two line^ for 
j = z + 2, . . . , 77, + z — 1 by considering Pn+r = Pr for r > 1. For every z G {1, . . . , n} 
in the case of a closed polygon and for every i E {1, . . . ,n — 1} in the case of an 
open polygon, we then obtain a value ji from the condition 

a] = min a]. (3) 

J' j=i+2,...,n+i-l ■' 



^The smaller one of the two complementary angles is taken in each case. 
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Once the point Pj^ has been selected by exploiting the illustrated criterion, we 
then establish the projective coordinate system {Xj, T^; Pj^} on the straight line Pj^Ti 
by calculating the projective representatives of Xj and Tj by solving 

jiXi + HiTi = Pj. 

for 7j and /ij. We obtain 

7i = Di^i/Di , jji = Di^2/Di , 

where 

A = det ( V Ai = det f P^-' V A2 = det ^ ^^'^ 



(4) 

/ ^ m G {0,1,2}. 

By Theorem [2] the point Ui is thus obtained as 

Ui = DiiXi — DioTi. 



3 Modification of the scheme for non-convex data 

If the input data are not convex, we segment them according to the following criteria 
in order to obtain piecewise convex segments. 

First, consecutive coUinear points are identified in the following way. 



By applying a dominant points selection algorithm like, for 
instance, the one in [20], we can easily detect the end points 
of a sequence containing at least 3 collinear vertices. Let us 
denote them by p° and p°. 



(5) 



Since in CAD applications as we have in mind linear features are usually intentional, 
we do not smooth the angle between a subpolygon consisting of (at least 3) collinear 
points and its neighbors. The insertion rule for these straight line subpolygons 
between the points and pf {k = 0, 1, 2,3,.. .) simply reads as: 



P2m= o(p' + Pti), ^=J, 1. (6) 



The remaining subpolygons do not anymore contain collinear segments. For these 
remaining subpolygons, inflection edges are identified by the following criterion, see 
also P[I6]. 

An edge p°p°_|_i is identified as inflection edge if the points , . 

p°_i and p°^2 lie in different half planes with respect to it. 

On an inflection edge we insert a new point, e.g., as midpoint of the edge corners, 
and we thus have a sequence of subpolygons without inflections. Next, we check 
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each of these subpolygons p° . . . p° for total convexity by the following criterion: 

If for every edge of the subpolygon all the points of the sub- 
polygon lie either on the edge or in the same half plane with 
respect to the edge, then the subpolygon is totally convex, 
otherwise it is only locally convex. 

If a subpolygon is only locally convex we divide it into two new subpolygons Pj • • • p° 
and p° . . . p°, where i = j + j. We repeat this procedure until we have a 

sequence of totally convex subpolygons, each of which we suppose to be composed 
of at least five points. 

Whereas in the interior of every convex subpolygon we apply the algorithm 
detailed in the previous section, we now define the method next to 

1. an inserted junction point on an inflection edge, which we refer to as inflection 
point, 

2. a junction point between two convex subpolygons, which we refer to as convex 
junction point. 

In the case of an inflection point p° let us denote the line of the inflection edge 
Pi'-iPiPi'+i by Cj, and the convex subpolygons meeting in p° by and s°j. We 
estimate a left and a right tangent in p?, 1°^ and 1°^ respectively, by applying the 
tangent estimation method from [1] to the five points p°_4, . . . , Pi_i, P^ of polygon 
s°i, respectively p°, . . . , p-+3, p-+4 of polygon s° .. 




Figure 5: Definition of an initial tangent 1^ in an inflection point p^. 



We then combine these two lines 1°^ and l^j for defining an initial tangent 1^ in p° 
(see Figure [5]): 

1° = + f^^r, , (9) 

where A° + /i° = 1, A°,/i° > 0. The pair (A°,yU°) thus plays the role of a shape 
parameter. The tangents in the other vertices of the polygons s^^ and s^^, and thus 
their new vertices, are calculated as described in the previous section by treating s^^ 
and s^ - separately. We obtain the new polygons s^^^^ and si -. Let us now describe 



9 



Figure 6: Tangent definition in an inflection point. 

how to obtain the tangents Igt^ in the point p° = Pgtj in the following iterations 
{k = 1, 2, 3, . . .), see Figure |6]for an illustration. 

Let 

St2''i = P2N-iP2^i and g^^2,. = p^.,p^.,+i (10) 
be the edges that are incident in p° = P2kp and 

lr,2H = ASr,2H^ 6^) and J^^,^ = ^{^l^H^ 6,) (11) 

their respective angles with the initial inflection edge e^. We then deflne the line ^^^^ 
by choosing that one of the lines gfgfej and g^2fei from (fTOj) yielding the maximum 
angle 

72% = max{7;;2.„7',2^J. (12) 
The tangent \^^^ in the point p° = Pg^^ is then deflned as 

where X'^^^ + /^2'=i ~ ■'■ and X^^k^, fi'2''i ^ ^- other vertices of sf^ and s^j we 

proceed as in section |2] for estimating the tangents; this allows us to calculate the 
new polygons s^f^ and s^+^ by separately applying the "convex" procedure from 
the previous section. 

In the case of a convex junction point p° we suppose our data to be such that 
the following condition hold^ (see Figure [7]): 

the intersection points Pi'_2Pi-i ^ Pi'Pi'+i and Pi„iP° H 

P°+iP°+2 lis in the same half plane with respect to the (14) 

line p°_]^p°^i as the point p°. 

Let us denote the convex subpolygons meeting in p^ by s^^ and s^^. As in the 
case of an inflection point we estimate a left and a right tangent in p°, and 

^This condition is guaranteed by a sufficiently dense sampling of the initial data points. 
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Figure 7: Situation in a convex junction point. 

l°j respectively, by applying the tangent estimation method from [1] to the five 
points p°_4, . . . , p°_i, p° of polygon sj' -, respectively p°, . . . , p°+3, p°+4 of polygon 
s° j. If the points p°_;^ and p°^]^ lie in different half planes with respect to 1°^ j 
respectively) we replace j respectively) by the line p°p°+i (Pi'-lPi respectively). 
We then combine these two lines and as in for defining an initial tangent 
1° in p°. The tangents in the other vertices of the polygons and s^.^, and thus 
their new vertices, are calculated as described in the previous section by treating 
and separately. We obtain the new polygons s^^^ and j. We then iterate 
this procedure and obtain the tangents l^k^ in the point p° = p^^k^ in the following 
iterations {k = 1,2,3,.. .) as 

where l^^^. and Ijjgfej is the respective left and right tangent in p^^^^ and A2fe^+/i2fej = 
'^2'=ii lA.^'i > (see Figure E] for an illustration). 




Figure 8: Tangent definition in a convex junction point. 
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In order to avoid a rather slow convergence of the newly inserted points towards 
the junction points (inflection points or convex junction points), see Figure [TOl we 
modify the point insertion rule for the points p^t+ii^i and P2fc+ij_,_i in the tangent 
triangles adjacent to the junction point p^^.. formed by the lines ^2''i-v ■^2'=i 
edge P2fei_iP2fci, and \^ki, ^^^i+i and the edge PafcjPsfcj+i respectively. The modified 
rule for the first and the last new point of a subpolygon reads as follows 

p = pt + am, p + a = l, p, cr>0, (16) 

where p stands for the new point Pgfe+ij.^ respectively Pgfe+ij+i; t designates the 
intersection point of the two tangents in the corner points of the corresponding edge 
and m is the mid-point of this edge, see Figure M 




Figure 9: Modified point insertion rule next to a junction point. 
This new "end point rule" avoids holes around the junction points (see Figure [TTI) . 



Pi / Pi * Pi 

/ * 

I /' 

I 

•I •/ •/ 

/ 

• It 4 • \ p 

-v.--.---* '-w.-..-.'* 

(a) fc = 1 (b) fc = 2 (c) fc = 3 

Figure 10: Example of slow convergence of the newly inserted points towards the 
inflection point pg. 
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® ® ^ 

•/ •/ 

•' 

-v.--.---® -w..--'* 

(a) fc = 1 (b) fc = 2 (c) fc = 3 

Figure 11: Application example of the new "end point rule" in equation f|T6l) to 
accelerate convergence of the newly inserted points towards the inflection point pg. 

4 The subdivision algorithm 

The above detailed procedure may be summarized in the following algorithm, which 
we apply to the given point sequence p° = (p° : i G Z). For the sake of clarity 
we will first describe the procedure to be used in the case of totally convex data 
(Subsection 14. ip . which constitutes an essential ingredient of the final algorithm for 
general non-convex data presented in Subsection 14.21 



4.1 Algorithm for totally-convex data 

Let p*^ = (pf : 2 = 1,... ,nfc) be the vertices of the totally convex polygon at the 
/c-level refinement. Hereafter we denote by the projective counterparts of the 
affine points pf, see section [21 

The algorithm that implements the function <y9 from ([T]) proceeds as follows in order 
to calculate the points p^;+\ = ^>{v\-2^ pf-i. Pf . P?+i. P?+2. v\+z \ P)- 

Step 1: Preprocessing 

a) For a closed polygon we set = P^, P^^ = P^.^, P^+i = Pf , P^+a = 
pfc 

-'2 • 

For an open polygon we set P^ = P^, P\ = Pf , P^+i = P4-4> ^4+2 = 
pfc 

For z = 1, . . . we then assign Qi = i^^_2, Q2 = Pj^Li, Q3 = Pt^ 

(^4 = Pj+i, = Pii+25 ^'^d apply formula for obtaining the tangent 
j^k pfc_ 

b) According to the angle criterion, we calculate the angles a^- for j = i + 
2, . . . , nfc + i — 1 and for i = 1, . . . , in the case of a closed polygon, 
while for i = 1, . . . , rifc — 1 in the case of an open polygon. We then select 
the value a*, satisfying condition ([3]). 
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Step 2: For a closed polyg on we set L^^_^_i — and we consider i — 1, . . . , ti/j, whereas 
for an open polygon we consider i = l,...,nfc — 1. We thus calculate the 
intersection points 

Tt = L\ A Ll, . 

Step 3: Calculate the lines (for i = l,...,nk for a closed polygon, and for i = 
1, . . . , TT-fc — 1 for an open polygon) 

as well as their intersection points 

Xf = ATf AA^ 



Step 4: Calculate the points P2i+\ as 



= DlX^ - Dl,T^ 



with D^^, according to 



4.2 Algorithm for non-convex data 

Step 1: Preprocessing 

We preprocess the data according to the criteria ([5]), ([7]) and ([8]) thus identify- 
ing subpolygons contained in a straight line and introducing and /or identifying 
inflection vertices and convex junction vertices within the initial data points 
of the remaining subpolygons yielding a sequence composed of totally convex 
subpolygons and "straight line" subpolygons. 

Step 2: • For a "straight line" subpolygon consisting of collinear segments we fix 
the tangents at its end points equal to the straight line passing through 
them and we apply everywhere the insertion rule ([6]). 

• For a totally convex subpolygon p° . . . = ... = p'^^.^^ . . . p^^.^ between 
the junction points p° and p° (where / > j + 4) we apply the algo- 
rithm of Subsection 14.11 in order to calculate the tangents in the points 
P2'=i+i' • • • ' P2''i-i points (with upper index k + 1) between 

P2fej+i ^^"^ P2'=i-i- According to the type of the subpolygon end points 
p'^kj and P2fep i.e., inflection point or convex junction point, we apply the 
rule (fT3l) for an inflection point and (fT5!) for a convex junction point in or- 
der to calculate the respective tangent, and apply rule (fT6l) for calculating 
the first respectively last new point, i.e., P2fe+ij+i respectively P2fc+i;_i- If 
p^kj or P2fc; is a junction point with a straight line subpolygon, we fix the 
tangent at such an end point equal to the straight line passing through it 
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and we define the new point in the adjacent triangle using the standard 
rule. Finally, if p^^.^ or p^^j coincides with an end point of the whole 
sequence, we proceed by computing the tangent at that location and the 
new point closest to it exactly following the same procedure we described 
in the open polygon case of Subsection 14.11 

5 Adaptive version of the scheme 

Subdivision curves are visualized by drawing a polyline on a level of refinement 
which evokes the impression of sufficient approximation of the given data. For a 
high quality rendering, the task is therefore to calculate and draw a level of subdivi- 
sion which is a visually sufficient approximation of the limit shape. The subdivision 
algorithm presented in Section H] provides a process of global refinement at every 
level. Therefore, when the starting polyline is highly non-uniform, the required 
level of refinement is determined by those locations which approximate the limit 
curve most unfavorably. Obviously, these may cause unnecessary fine subdivisions 
at other locations of the curve, thus leading to an unreasonable resource demanding 
algorithm. 

To overcome this problem we propose an adaptive version of the subdivision algo- 
rithm previously described. Adaptive subdivision is achieved by applying the mech- 
anism of subdivision only locally, i.e. only at those locations of the initial polyline 
that are not approximated with the desired quality. The decision where high res- 
olution refinement is needed, strongly depends on the underlying application. As 
concerns our algorithm, adaptivity may be controlled either by the user or by an 
automatic criterion. In fact, the user may specify which portions of the polygon 
should be subdivided or the process may be automated by controlling whether the 
length of an edge is greater or not than a specified threshold. Only in the positive 
case we insert a new point in correspondence of the considered edge. Of course, 
besides improving the visual quality of the limit curve, the adaptive version of the 
scheme reduces the computational cost of the algorithm. 

In the remainder of this section we take highly non-uniform polylines and compare 
our adaptive refinement algorithm with the basic one. Figure [12] shows the com- 
parison between the refined polylines obtained by applying the two algorithms to 
a totally convex set of points, while Figure [M] compares the two algorithms on a 
highly non-uniform polyline with convex junction points and infiection points. 

Figures [T5] and [T^ illustrate the corresponding discrete curvature plots. It can be 
easily seen that the use of the adaptive refinement scheme results in a considerable 
improvement of the curvature behaviour. 
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(a) 



(b) 



(c) 



Figure 12: Application example of the adaptive algorithm: (a) starting polyline; (b) 
refined polyline after 6 steps of the basic algorithm; (c) refined polyline after 6 steps 
of the adaptive algorithm. 



(a) (b) 

Figure 13: Comparison between the discrete curvature plots of the refined polylines 
in Figure [H (b) and [H (c). 




(a) (b) (c) 

Figure 14: Application example of the adaptive algorithm: (a) starting polyline; 
(b) refined polyline after 6 steps of the basic algorithm; (c) refined polyline after 6 
steps of the adaptive algorithm. The labels in Subfigure (a) denote the end points 
of consecutive subpolygons: while pj is a convex junction point, and p'j's 
infiection points. 
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Figure 15: Comparison between the discrete curvature plots of the refined polyhnes 
in Figure [H (b) and [H (c). 

6 Properties of the scheme and smoothness anal- 
ysis 

By construction the presented scheme enjoys the properties of being shape preserving 
and conic reproducing as summarized in the following proposition. 

Proposition 1. The presented interpolatory curve subdivision scheme 

a) is shape preserving, i.e., if the starting point sequence p° = (p° : i = 
1, . . . , uq) consists of convex, straight line and concave segments, then all gen- 
erated subsequent point sequences p'^ = (p^ : i = 1, . . . , rik), for /c = 1, 2, 3, . . . 
respect the same behavior; 

b) reproduces conic sections, i.e., if the starting point sequence p° is sampled from 
a conic section c, then the limit curve coincides with c. 

Proof. a) In the preprocessing step the initial point sequence is segmented into 
straight line segments and totally convex segments, which are not changed 
during the subdivision procedure. A straight line subpolygon is reproduced 
as such. For a totally convex subpolygon p'^ = (pf : i = jk,---,h) the 
algorithm first generates a line If in every point pf. By construction (see 
Section [2] and (fT3l) . (fT5|) ) these lines form a convex delimiting polygon for the 
points of the sequence p'^"'"^ of the next subdivision level, see Figure [2l The 
fact that by construction the point P2i'+i is contained in the triangle formed 
by the points p*^, tf , p*^_;^ (P2i+i ^ ^(Pi ^ '^i ^ Pi+i)) guarantees convexity of the 
sequence p'^"'"^. 

b) If the point sequence p'^ comes from a conic section c, then the lines If gen- 
erated in the preprocessing step are the tangents of c in the points pf. By 
Theorem [2] the constructed points P2i+i he on c. □ 
The special set-up of the presented subdivision scheme allows to obtain a result 
on the smoothness of the limit curve as formulated in the following proposition for 
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which we suppose our data not to contain any subpolygons that consist of consecu- 
tive colhnear segments. The transitions from and to such subpolygons are purposely 
by construction, and these "straight line" subpolygons are exactly reproduced 
by the subdivision algorithm. We can thus restrict our attention to polygons not 
containing any straight line subpolygons. 

Proposition 2. a) For the presented subdivision scheme the polygon series p'^ 
converges to a continuous curve. 

b) The limit curve of the presented subdivision scheme is of continuity class G^. 

Proof. The proof is formulated for the non-adaptive version of the subdivision al- 
gorithm where every polygon edge is replaced by two new edges in every step. For 
the sake of simplicity we omit the details for the adaptive case, where the indices 
change but not the general idea. We proceed by demonstrating the C° and con- 
tinuity in two steps. First, we consider (locally) convex segments (including convex 
junction points), and then we treat inflection points. The locally convex segments 
are composed of several totally convex segments joined together at convex junction 
points. 

We introduce the following notation. See Figure [I6] for an illustration. 




Figure 16: Notation used throughout the proof. 



intersection point of the edges pf_ipf and p^^xP^^_2 

inner angle in pf of the triangle A(pf^_^pfqf) 
/3f inner angle in pf of the triangle A(p^^_^pftf) 
vrf inner angle in pf of the triangle A(p^^_;^pf P2j^\) 
h'l height of the triangle A(pf_^^pftf ) from onto the edge pfpf+i 

height of the triangle A(pJY;^pf P2jt^\) from P2t+i onto the edge pfpj^+i 

length of the edge pf pj+i 

length of the line segment pf 
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C° continuity for (locally) convex segments: 



In order to show the continuity of the limit curve we compute the distance d!' between 
the polygon p'^"'"^ and the polygon p*^. According to the above notation we define 
d!^ = maxj{(i^}. By constructioijfl the point t*^ lies inside the triangle A(p^^^pfq*^) 
and the point ^2t+i ^i^^ inside the triangle A(p*Y]^p^tf ). Thus (see Figures [T6 | [T7|) 



< vrf < < 



(17) 



and 



P2i ^ Pi 




k _ k+l 
Pi+1 — P21+2 



Figure 17: Relation between the angles pf and P2 



k+l 



Let pj" be a point that appears in the ko-th iteration. Since the presented scheme 
is interpolatory we have 



„fco _ „fco+l _ „fco+2 
From (fTSl) we obtain in the point p^°: 



(19) 



and by fjTTj) we have 



(20) 



■^In a convex junction point this is guaranteed by condition (|14p and the respective tangent 
definition. 
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where the e/s are constants such that < < 1. Let e := maxj{ej} and 
maxjlpy"}, where < e < 1. Then, 



P2"j 


< 


^npko 


oko+n 
P2"j 


< 


^npko 


fco+n 
^2"i 


< 


^npko 



(21) 
(22) 
(23) 



In the triangle A{p^_^_lP^p2t+l) have 



sinfTT,^) 



,fc+i 



C2i 



and thus 

4°+" = 4°+7Sin(7r2'«+"). (24) 

Since the angles become smaller in every step and the longest edge of a triangle 
is the one opposite to its biggest inner angle there exists an index k' such that for 
all k > k' the circle centered in with radius contains the triangle A(p^_^]^p^q^'). 
Thus there exists an index Uq such that for all n > Uq and for all j we have 



c, 



,A;o+"-+l ^ „fco+"0 + l ^ mavl n^O+no + lX _. „fco+no + l 



Since sin(x) < x for x > 0, we deduce from (!24|) and (!23|) for n > Uq that: 

jko+n ^ fco+no+l„fco+" 



Since this holds for all j we have 

^ ^fco+no+lgfc-fcopfco ^ (25) 

where A; > /cq + no. The polygons {p*^} thus form a Cauchy sequence and this se- 
quence of polygons converges uniformly. Each polygon being piecewise linear, the 
limit curve is continuous. 



continuity for (locally) convex segments: 

Let 

h'' = max{hi} . 

i 

For the triangle A(pf_^^pft*^) it holds 

sin(/3f ) = f . 

By the analogous reasoning as above for we obtain for (since there exists an 
index k' such that for all k > k': /f < c^): 

Lfc ^ „ko+n k-ko ko 
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where k > ko + n for a certain value of n and 

c'=o+^ = max{4f } . 

In every iteration level the presented subdivision scheme constructs the new points 
by sampling them from a -continuous conic spline T''. By its Bezier construction 
the distance of the conic segment of T'' corresponding to the edge pf+i is bounded 
by hk- The sequence of conic splines {T''} thus converges to the same limit curve as 
the sequence of polygons {p'^}. 

It thus remains to be shown that the sequence of tangents {l'^}, where l'^' = (If : z G 
Z) converges uniformly. Let ^2i^^ be the angle between the tangents If and Igj^^ in 
pf. Then, 

see Figure fT8l 




Figure 18: Illustration of the uniform convergence of the tangent sequence. 
Thus by ([22]) and (ES]) with i = 2"j and k = k^ + n: 

l^-fco+n+ll ^ \nko+n+l\ . \ ko+n\ . \nko+n\ 

< e"(e + 2)|p'=o|. 

The set of tangents {l'^} thus forms a Cauchy sequence and this sequence converges 
uniformly. The limit curve is therefore of continuity class G^. 

C° and continuity for inflection points: 
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In an inflection point p° = Pg^^ we deflne 



k 


— ^(l2fe-ij'S2fcj 








= Z(l2fc_i^, Igfcj) 


^^2^1 


= ^(l2fej'Sz^2fci) 


Pr,2*=i 


~ ^{^2''ii Sr,2'=i) 



where Igj-, ggfej, gf2'=j' ^^"^ Sr2'=i deflned as in section El see Figure fT9l We have 

4. = 4-. - , (26) 

with 72fc. from (fT2|) . 

\\4 















^^^.^P2'',+l 
















\ 1* 





Figure 19: Illustration of C° and continuity for inflection points. 
By construction these angles satisfy the following inequalities (see Figure [19]): 



o<e2% < 4^, (27) 

0<4, < 4-., (28) 

0<72% < 7ji\.' (29) 

0</3^2+i,, < ^3l,2''^^ (30) 

0</3,%V., < (31) 



The relations (130|) and (13T!) imply that the tangent triangles adjacent to the in- 
flection point p°, used for the determination of the new points next to p°, become 
continuously flatter, thus yielding continuity in p°. 
Furthermore, we have 

^2H = ^2>^-H - I2H > ^2H - 72fe+i, = t^2fe+i* > • (32) 
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By introducing an e with < e < 1 we thus have 



and consequently 

The angle between two consecutive tangents in an inflection point p° thus continu- 
ously decreases yielding a continuous inflection joint in the limit. □ 



7 Numerical examples 

In this section we want to illustrate the performance of the basic and adaptive 
subdivision algorithms presented in Sections H] and respectively. Generally, when 
the starting points are uniformly spaced we use the flrst proposal, while in the case 
of irregularly distributed vertices we apply the second one. 

As concerns the forthcoming examples, we start by applying the subdivision al- 
gorithm of Subsection 14.11 to totally convex closed and open polylines with nearly 
uniform edges (Figures [2U | 1^2]) . and successively we exploit the adaptive version of 
the scheme in the case of polylines with highly non-uniform edges (Figures |2I1 125]) . 
As it appears, the generated curves are always convex and visually pleasing. 
Then, concerning Figures [2l] and [251 the goal is to illustrate the conic precision 
property of the proposed algorithms in both the uniform and non-uniform cases. 
Moving from top to bottom, the four plots in the flgures have been generated by 
repeated application of the subdivision scheme to points sampled from a circle, an 
ellipse, a parabola and a hyperbola. 

We then continue by showing that the curves computed through the algorithms pre- 
sented in Subsection 14.11 and Section [5] are really artifact free. In fact, although a 
limit curve can be apparently artifact free, it is hard to tell from the display on the 
screen if it is acceptable or not. Two curves may look very similar on the screen, but 
their curvature plots may reveal important differences. The most commonly used 
tool for revealing signiflcant shape differences is provided by the curvature comb of 
the curve. In pictures [261 127] we have used the graphs of the discrete curvature combs 
of the reflned polylines to show that the limit curves generated by our algorithms 
are indeed artifact free. 

In particular, if we compare the results we get by reflning the polyline in Figure 
[271(a) for data, that do not come from a conic section, through our adaptive algo- 
rithm and through the subdivision algorithms in [6], [H] and [25] (Figures [27[(b). 
(c) and (d) respectively), they are only apparently very similar. Yet their curvature 
combs reveal substantial differences showing that neither every non-linear nor every 
non- uniform subdivision scheme is indeed artifact free (see Figures [57] (e)-(f)-(g)- 
(h)). 

We close this section by illustrating the results of the subdivision algorithm of Sub- 
section 14.21 In the flrst example we take the D-shape polyline of [9] in order to 
show the ability of the scheme to reproduce collinear vertices (see Figure [251) . In the 
second and third examples we apply the generalized subdivision scheme to a closed. 
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respectively open, sequence of non-convex data to illustrate its continuity and 
shape-preserving interpolation properties (Figures 129 } |30|) . The last three examples 
in Figures [3T] and |32] deal with more complex polylines that respectively represent 
the cover of a mobile phone, Micky Mouse face and a bottle opener. 
The shapes in Figure El] are designed by a collection of rectilinear and convex seg- 
ments where the most of the latter ones have been sampled from conic sections. 
The shape in Figure [32] is made of 4 independent closed polygons, some of which 
are non-convex. The data of the first example in Figure [31] and that of Figure [32] 
are courtesy of the CAD company think3 (www.think3.com). In all the considered 
experiments the proposed subdivision scheme turns out to work very well and clearly 
manifests all its characteristic features described in the previous sections. 




Figure 20: Application examples of the subdivision algorithm of Subsection 14.11 to 
nearly uniform closed polylines. From left to right: points at 1st and 2nd level of 
refinement; refined polyline after 6 steps of the algorithm. 




Figure 21: Application examples of the adaptive subdivision algorithm of Section [S] 
to highly non-uniform closed polylines. From left to right: points at 1st and 2nd 
level of refinement; refined polyline after 6 steps of the algorithm. 
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Figure 22: Application examples of the subdivision algorithm of Subsection 14.11 to 
nearly uniform open polylines. From left to right: points at 1st and 2nd level of 
refinement; refined polyline after 6 steps of the algorithm. 



i 

Figure 23: Application examples of the adaptive subdivision algorithm of Section 
[5] to highly non-uniform open polylines. From left to right: points at 1st and 2nd 
level of refinement; refined polyline after 6 steps of the algorithm. 

8 Conclusions 

Even though in the last years important steps forward have been taken both in the 
construction and analysis of interpolatory subdivision schemes |2l], several problems 
are still open and need to be tackled in order to increase the strength and popularity 
of subdivision in more and more fields of application. 

First of all, unlike the non-interpolatory subdivision schemes, the interpolatory 
ones usually generate shapes of inferior quality because, if applied to points with 
an irregular distribution, they provide a limit curve with more convexity changes 
than the starting polygon. Since, in several applications it is important to guarantee 
shape preservation, in this paper we have described a new interpolating subdivision 
algorithm enjoying this important property. 

Because in CAGD it is also often necessary to have schemes able to generate clas- 
sical geometric shapes, we have enriched our interpolating subdivision scheme with 
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Figure 24: Uniform case: reproduction of conic sections from uniform samples by 
applying the subdivision algorithm of Subsection 14.11 The dotted blue line is the 
conic section to be reconstructed. From left to right: points at 1st and 2nd level of 
refinement; refined polyline after 6 steps of the algorithm. 

the capability of including the exact representation of all conic sections. The differ- 
ent methods of the literature [U [5l [22] give solutions only when the assigned points 
have a regular distribution. These linear subdivision schemes use non-stationary 
refinement rules associated with functional spaces defined via the union between 
polynomial and exponential functions with a free parameter. 

The idea we have explored in this paper is to provide a subdivision scheme in 
which the conic section reproduction is obtained by adaptive geometric constructions 
on the given points. The advantage of doing this is that the presented non-linear 
scheme is able to adapt itself to any data configuration, i.e., to arbitrary irregularly 
distributed point sequences. Due to the underlying construction, the properties of 
shape preservation, conic reproduction as well as the proof of continuity of the 
limit curve, follow straightforwardly. The method has been illustrated by several 
significant examples. 
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Figure 25: Non- uniform case: reproduction of conic sections from non-equispaced 
samples by applying the adaptive subdivision algorithm of Section [51 The dotted 
blue line is the conic section to be reconstructed. From left to right: points at 1st 
and 2nd level of refinement; refined polyline after 6 steps of the algorithm. 
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(a) (b) (c) (d) 




(e) (f) (g) (h) 

Figure 27: Comparison with the hnear, non-uniform scheme in [6] and with the 
non-hnear schemes in [H] and [25]. First row: refined polyhnes obtained after 6 
steps of (a) our adaptive algorithm of Section [5l (b) algorithm in [6] with chord 
length parameterization; (c) algorithm in [H] with chord length parameterization; 
(d) algorithm in |25]. Second row: corresponding curvature combs. 




Figure 28: Application example of the subdivision algorithm of Subsection 14.21 to a 
closed sequence of data containing coUinear vertices. The labels in the figure denote 
the end points of consecutive subpolygons. From left to right: points at 1st and 2nd 
level of refinement; refined polyline after 6 steps of the algorithm. 
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[17] Kuijt, F., 1998. Convexity preserving interpolation - stationary nonlinear sub- 
division and splines. PhD Thesis, University of Twenty. 

[18] Kuijt, F., Van Damme, R., 2002. Shape-preserving interpolatory subdivision 
schemes for nonuniform data. J. Approx. Theory 114, 1-32. 

[19] Marinov, M., Dyn, N., Levin, D., 2005. Geometrically controlled 4-point in- 
terpolatory schemes. In: Dodgson, N.A., Floater, M.S., Sabin, M.A. (Eds.), 
Advances in Multiresolution for Geometric Modelling, Springer- Verlag, pp. 301- 
315. 

[20] Marji, M., Siy, P., 2003. A new algorithm for dominant points detection and 
polygonization of digital curves. Pattern Recognition 36, 2239-2251. 

[21] Pascal, E., 1910. Repertorium der Hoheren Mathematik. Teubner, B.G., 
Leipzig/Berlin. 

[22] Romani, L., 2009. From approximating subdivision schemes for exponential 
splines to high-performance interpolating algorithms. J. Comput. Appl. Math. 
224(1), 383-396. 



30 



V.' : 1 
ii i 

1 ' ' 


1 1 






i i 








1 •-■•■^ •'•-^ 1 
•^•-^ •^**^ 

•- — ••- — • •- — • ' 

' ' 

♦ - — '^M 

1 1 


ooo 
ooo 
ooo 
ooo 




0000 
0000 
0000 


1 i 










Figure 31: Application examples of the subdivision algorithm of Subsection 14.21 to 
complex data. Left: starting polylines. Center and Right: refined polylines obtained 
after 7 steps of our algorithm. Data of first row: courtesy of think3. 




Figure 32: Application example of the subdivision algorithm of Subsection 14.21 to 
complex data. Left: starting polylines. Center and Right: refined polylines obtained 
after 7 steps of our algorithm. Data are courtesy of think3. 
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